Truncation of periodic image interactions for confined systems 
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First principles methods based on periodic boundary conditions are used extensively by materials 
theorists. However, applying these methods to systems with confined electronic states entails the use 
of large unit cells in order to avoid artificial image interactions. We present a general approach for 
truncating the Coulomb interaction that removes image effects directly and leads to well converged 
results for modest-sized periodic cells. As an illustration, we find the lowest-energy quasiparticle 
and exciton states in two-dimensional hexagonal GaN sheets. These sheets have been proposed as 
parent materials for single-walled GaN nanotubes which may be of interest for optoelectronics. 
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First principles methods, especially those based on 
Density Functional Theory (DFT) 0, |2| , are used widely 
to calculate physical and chemical properties. One very 
popular approach uses plane wave (Fourier) basis sets, 
periodic boundary conditions, and pseudopotentials 
Plane waves provide a general, orthonormal basis with 
controllable convergence. They also are natural for phys- 
ical situations where electronic states extend in all di- 
rections (e.g. bulk materials). Mature, well optimized 
and parallelized plane-wave software packages are pub- 
licly available ;4]. In addition to electronic ground-state 
properties from DFT, quasiparticle and optical excita- 
tions can be calculated with methods that build on this 
theoretical and computational edifice |E IE S B IS ^| ■ 

Many systems of interest have spatially confined elec- 
tronic states. Molecules are classic cases, but nanostruc- 
tures such as nanotubes, nanowires, or sheet-like sys- 
tems are of current interest. Often, it is the confine- 
ment of electronic states that leads to an unusual be- 
havior. Unfortunately, studying these structures with 
plane waves requires the use of large unit cells with re- 
gions of empty space surrounding the structures. This 
is necessary to ensure that periodic copies do not inter- 
act: one usually aims to study the isolated structure and 
not the artificial array. The electronic wave functions 
are not a problem because they decay exponentially into 
the vacuum regions and converge rapidly with increasing 
periodic separation. Rather, the main difficulty stems 
from the long-range Coulomb interaction. It can lead 
to long-range monopolar or dipolar image interactions 
that fall off slowly (algebraically) with periodic separa- 
tion 0, ^|. In practice, one has to model the image 
interactions and extrapolate to infinite separation j^j . 

Ideally, one should eliminate Coulomb interactions 
among periodic copies from the start. Truncation of 
the Coulomb interaction was used for molecular sys- 
tems and carbon nanotubes [T^ IT^ . Here we provide 
a general scheme for truncating the Coulomb interac- 
tion in systems having both periodic and confined spa- 
tial directions. The scheme should interest researchers 
investigating confined systems with plane wave meth- 
ods or, more generally, periodic boundary conditions be- 
cause it leads to rapid convergence versus periodic sepa- 
ration. Our illustrative calculations focus on calculating 
electronic excitations where image effects are most pro- 



nounced However, the approach is equally useful 

and easy to implement for ground-state methods. 

We begin with ground-state calculations using plane 
waves. The component of the total energy incorporating 
long-range Coulomb interactions is the Hartree energy, 



Eh = d?r I d\' p{r)v,{r - r')p{r') 



i^|p(G)p5,(G), 



where Vc{r) is the Coulomb kernel Vc{r) = l/|r| for 
the untruncated case, and p{r) is the total charge den- 
sity. The G are the reciprocal lattice vectors of the pe- 
riodic cell, and /5(G) is the Fourier transform of p{r). 
The Fourier transform of the untruncated interaction is 
Vc{k) = 47r/fc^. (The G = term is excluded due to 
a compensating background 0.) The remaining terms 
in total energy, the kinetic and exchange-correlation en- 
ergies, are short-ranged for the standard semi-local ap- 
proximations such as the LSDA or GGA We want 
to replace Vc{r) by a function that explicitly zeroes in- 
teractions between periodic images. Obviously, the re- 
placement must yield a Fourier transform Vc{k) that is 
well-defined for fc ^ 0. 

For an isolated atom or molecule, we truncate in all 
three spatial directions. A simple spherical truncation 
with radius Tc was used previously |l2lll4| . 
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The interaction v^^{k) is finite for any k. The issue is 
to choose Tc- One takes Tc to be larger than the "size" 
of the molecule d (the region where p(r) is appreciable). 
Thus intramolecular interactions are calculated correctly. 
Simultaneously, rc must be smaller than the shortest dis- 
tance L between periodic images minus d to prev ent ' 
odic image interactions, leading to d < < L—d \ 
A practical solution is to set = L/2 and increase L to 
achieve convergence. 

The smoothness of Vc{k) is crucial for calculations of 
electronic excitations. Green's function-based methods 
describing quasiparticle and excitonic dynamics require 
the screened Coulomb interaction matrix WcciQ) where 
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g is a wave vector in the first Brillouin zone 0, . (We 
will truncate interactions in the confined directions so q 
only has components along the non-confined, i.e. physi- 
cally periodic, directions.) W is calculated through the 
dielectric matrix e. 
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GG' 



WGG'{q) = e 
f^GG'iq) = Sgg 



-G') 

G) XGG'iq) , 



where x is the irreducible polarizability typically evalu- 
ated within the random-phase approximation (RPA) 
(The frequency dependences of W, e, and x ^re sup- 
pressed for clarity.) Calculations of excited state proper- 
ties require sums over the entries of W and an integral 
over q. With spherical truncation, vl^{k) is well behaved 
making for straightforward calculation of W . For the 
usual untruncated case, Vc{k) = 47r/fc^ diverges as fc ^ 0, 
but the key divergence in W is confined to the "head" 
element G = G" = when g ^ H [H [H. The cal- 
culation and inversion of e is handled separately for this 
special entry .5]. In brief, divergent behavior in W is 
confined to one known entry that is treated "by hand" . 

We derive analogous truncations for the two remaining 
cases: (1) sheet-like geometries with one confined and two 
periodic directions (e.g. graphene or thin films), and (2) 
wire-like geometries with one periodic and two confined 
directions (e.g. nanotubes or nanowires). 

We begin with the sheet geometry. The confined di- 
rection is along the z axis, and the system is periodic in 
the xy plane. We choose a truncation length and set 
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We require d < Zc < Lz where d is the the extension of 
electronic states into the inter-sheet vacuum region and 
Lz is the periodicity along z. The Fourier transform is 
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siiv{kzZc) - cos(fc^Zc) 



where k^y — {k'^ + kyY^^. As a check, in the short- 
wavelength limit kxyZc ^ 1, we recover the untruncated 
interaction. For arbitrary Zc, this formula diverges when 



for any nonzero kz and is clearly unsuitable. 



However, the choice Zc = Lzl2\s special. For the sheet 
geometry, reciprocal lattice vectors have kz — 2TTnz/Lz 
for integer Uz- Thus sm(kzZc) = for all Uz when Zc = 
Lzl'2.. We arrive at the simplification 



f(fc) 
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The Coulomb interaction is now finite for any k^y when 
kz ^ ^- For kz = 0, it diverges as A-'nZclkxy as k^y — > 0. 
As expected, the divergence is milder than the untrun- 
cated 3D case. Our choice of Zc has confined the di- 
vergence to the single wave vector fc = 0. This simple 
replacement for Vcik) allows for straightforward ground- 
state calculations. (The minor required changes for ex- 
cited states are outlined below.) 



We proceed to the wire geometry. The periodic direc- 
tion (wire axis) is along z, and the system is confined in 
the xy plane. We seek a truncation of the form 



9ix,y) 



(4) 



The shorthand 9{x, y) returns zero unless x and y lie in a 
finite region surrounding the wire that is specified below. 
The Fourier transform along z gives 

vrik)^ fdxfdy0{x,y)2Ko{\kz\p)cos{k,x+kyy). (5) 



Here, Kq{z) is the modified Bessel function that diverges 
as -ln(z/2) for z and p = {x'^ + y^)^/^. The xy 
integral is of finite extent so divergences in v^^{k) orig- 
inate from Kq as kz — > 0. Using the asymptotic form 
Ko{z) = — Inz + 0(z'''), we isolate the divergent term 



-2 ln(|A;2 1) J dx J dy 9{x, y) cos{kxX -I- kyy) . 

For any wave vector with kxy ^ 0, we make this di- 
vergence vanish by setting 0(x^ y) 7^ in exactly one 
Wigner-Seitz cell centered on the wire in the xy plane. 
The above xy integral is then zero when k^y ^ because 
the projection of k in the xy plane is a reciprocal lattice 
vector. This make all entries of vf^{k) finite when fc 7^ 0, 
and only fc = diverges. For an arbitrary unit cell in the 
xy plane, no further simplification is possible and w^*(fc) 
must be computed numerically. 

This choice is superior to what one may have originally 
used: cylindrical confinement 



qcyl 
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This cylindrical choice gives 
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2 {l + k^yPcJi{k^yPc)KQ{\kz\pc)- 



Although we have a closed form, for any kxy ^ the Kq 
term diverges as kz 0. (This is because a 2D unit cell 
in the xy plane can not have a cylindrical cross section.) 
Many entries Wgg' (?) are divergent, and this leads to nu- 
merical difficulties. In practice, we integrate the entries 
WcG'iq) over g, and the integrals are finite. However, 
performing the integrals requires very dense q sampling 
In contrast, our mathematically justified choice con- 
fines the divergent behavior to the single element fc = 
in {)™*(fc). As before, this leads to one divergent entry 
WoQ{q) that is handled separately. 

We illustrate our findings by computing the quasipar- 
ticle and optical band gaps of two-dimensional hexagonal 
sheets of GaN. These sheets have been investigated with 
but not with more accurate excited- 
state methods. In analogy with graphene and carbon 
nanotubes, these GaN sheets may be the parent materi- 
als for single- walled GaN nanotubes, systems of possible 
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FIG. 1: Static inverse dielectric constant eggiq) versus mo- 
mentum transfer q for GaN sheets. The periodic separa- 
tion between sheets is 14 a.u. and Coulomb truncation was 
used. Red diamonds are the first principles results using 
the RPA. The solid line is a fit of the form l/eQQ{q) — 
1 -1- 'y\qfvc'^{q) exp(-Q|g|) with 7 = 2.2 and a = 2.7 a.u. 



interest for their optical and luminescent properties. The 
lack of accurate excited-state theory is a serious short- 
coming because DFT is known to fail in predicting these 
very properties |^ . 

We describe the ground state using DFT-LDA with 
nonlocal norm-conserving pseudopotentials and a plane- 
wave cutoff of 35 Ryd 3, 22] . The hexagonal sheets have 
a primitive lattice vector of 5.92 a.u. in the xy plane. 
The supercell is created through k-point sampling of the 
primitive cell with a regular Monkhorst-Pack grid [23l |. 
Within DFT-LDA, the system is insulating with an in- 
direct band gap of 2.4 eV between F and K. We concen- 
trate on the lowest-energy direct transition which is at F 
and has an energy of 3.1 eV within DFT-LDA. To calcu- 
late the excited states, we employ the GW approxima- 
tion to the self-energy for quasiparticles and the Bethe- 
Salpeter Equation (BSE) for optical properties 0, H^ . 
The dielectric matrix has a cutoff of 20 Ryd. Sums over 
unoccupied states in the GW calculations extends to 22 
Ryd. A static electron-hole interaction is used in the 
BSE. All remaining parameters are chosen to converge 
excitation energies to better than 0.05 eV. 

We first consider the quasiparticle band gap at F. Ta- 
ble ^ shows the GW correction to DFT-LDA as a func- 
tion of k-point sampling and the separation of periodic 
copies along z when no Coulomb truncation is employed. 
The idea is to densely sample in the xy plane to sim- 
ulate an infinite sheet, consider different periodic sepa- 
rations, and try to find the limit of infinitely-separated 
sheets. This corresponds to n x n x 1 k-point sampling 
of the primitive cell. Unfortunately, convergence with re- 
spect to k-sampling is slow. Even for very large supercells 
(14 X 14 X 1), the results are not converged to the typical 
required accuracy of 0.1 eV. This difficulty is not related 



No truncation: GW gap correction (eV) 



k-point 


Sheet separation Lz (a.u.) 


sampling 


10 


14 


18 


4x4x1 


1.57 


1.80 


2.01 


8x8x1 


1.98 


2.05 


2.10 


12 X 12 X 1 


2.21 


2.24 


2.24 


14 X 14 X 1 




2.33 


2.32 






8x8x1 


1.98 


2.05 




8x8x2 


1.78 


1.99 




8x8x3 


2.06 


2.06 





TABLE L GW correction to the DFT-LDA minimum direct 
band gap at F for GaN sheets as a function of k-point sam- 
pling and periodic sheet separation. No Coulomb truncation 
is employed. See text for details. 



With truncation: GW gap correction (eV) 



k-point 


Sheet separation Lz (a.u.) 


sampling 


10 


14 


18 


4x4x1 


2.39 


2.60 


2.83 


8x8x1 


2.53 


2.55 


2.58 


12 X 12 X 1 


2.53 


2.53 


2.52 


14 X 14 X 1 




2.53 


2.52 



TABLE II: Same as Tabled] but with Coulomb truncations. 



With truncation: Exciton binding energy (eV) 



k-point 


Sheet separation Lz (a.u.) 


sampling 


10 


14 


4x4x1 


2.14 


2.11 


8x8x1 


1.46 


1.40 


10 X 10 X 1 


1.37 


1.32 


12 X 12 x 1 


1.33 


1.30 


14 X 14 X 1 




1.30 



TABLE III: Binding energy for the lowest optically active 
exciton in the GaN sheets. Coulomb truncation is employed. 



to the k-sanipling along z: as the Table shows, k-point 
sampling along z (n x n x m sampling) does not improve 
the situation. Fundamentally, the long-range image in- 
teractions are hindering the convergence. 

By contrast, our truncation method yields well- 
converged results with modest effort. Before present- 
ing the results, we summarize the minor changes re- 
quired to treat the diverging element Woo (9) in the GW 
and BSE calculations. The plan is identical to the case 
with no truncation: we replace the divergent element 
by an appropriate average over q. In the GW calcu- 
lations, we average Woo (9) over the first Brillouin zone 
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of the supercell @ . The average can be performed eas- 
ily by using a model fitted to the ab initio results. We 
have Woo (9) = ^mil)^c^(l) where we have calculated 
^00^(9) ^ S''i'i °f nonzero q vectors. In addition, 
simple analytical considerations show that for insulators 
l/eQQ{q) = l+v^'^{q)f{q) where f{q) is smooth and has a 
Taylor series beginning at quadratic order "iG] . A simple 
form f(q) = j\q\'^ exp{—a\q\) with two positive parame- 
ters a and 7 is sufficient. See Figure^for an illustration. 
The fit need be good only in the region of small q since 
the Brillouin zone of the supercell is a small region close 
to q = 0. For the BSE, similar types of averages are 
needed 0|, and we also use the model function. 

Table m shows that our Coulomb truncation method 
leads to rapid convergence. We conclude that the GW 
correction for an isolated sheet is 2.5 eV, making for a 
total direct quasiparticle band gap of 5.6 eV at F. Ta- 
ble mil displays the convergence of the binding energy 
of the lowest-energy exciton calculated with the BSE. 
This exciton is optically allowed. We conclude that the 
binding energy is 1.3 eV. Thus the lowest exciton for 



the sheet is has an excitation energy of 4.3 eV and is 
optically bright. However, some caution is required be- 
fore extrapolating this result to the optical properties of 
GaN nanotubes. In semiconducting carbon nanotubes, 
dark (diople-forbidden) excitons are the lowest in energy 
likely leading to detrimental reduction of luminescence 
efficiency 24] . Our preliminary calculations on GaN nan- 
otubes point to the same problem. 

In summary, we derive a scheme for truncating the 
Coulomb interaction in plane-wave calculations of sys- 
tems with confined electronic states. The method 
excludes periodic image interactions and yields well- 
converged results for modestly sized periodic cells. We 
illustrate the method by using it to calculate the lowest 
energy optical excitations in two-dimensional hexagonal 
GaN sheets. 

Computer resources were provided by Yale High Per- 
formance Computing. This work was facilitated by the 
U.S. Department of Energy's Computational Materials 
Science Network (CMSN). 
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